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A lecture  presented  at  Instltut  de  Recherche  d' Informatique  et  d'Automatique 
in  Rocquencourt , France,  on  May  11,  1978.  This  lecture  is  included  In  the 
publication  Semlnaires  IRIA  Analyse  et  Contrflle  de  Systemes  - 1978. 


1.  Introduction 

The  need  for  the  reformulation  of  hydrodynamics  presented  here  arose  in  con- 
nection with  the  study  of  the  hydrodynamic  free  surface  problem.  As  the  problem 
is  conventionally  formulated,  one  looks  for  a solution  (u(x,t) ,P(x,t) ,(^_(t))  of 
the  initial  value  problem  for  the  equations 


u + u • Vu  ■ - ^ VP-gk  , xe(/^(t)  , 0 < t < T, 

p0 

V • u - 0 , xe  (t) , 0 < t < T, 


(1.1a) 

(1.1b) 


where  u is  the  velocity  field  and  P is  the  pressure.  pQ  is  the  density  of 
water,  g is  the  gravitational  constant,  and  k is  a unit  vector  in  the  z-direction. 


The  initial  conditions  for  (1.1)  are 

u(x,0)  - uQ,  (1.2a) 

$(0)  -Q0-  (1.2b) 

Boundary  conditions  are 

P ■*  - PQgz  as  z -*■  - « , (1.3a) 

u • n » • n on  3(^g(t)  , (1.3b) 

P - 0 on  3$f(t)  , (1.3c) 


Separate  copies  of  this  paper  are  being  sent  to  addresses  on  the  distribution 
list  for  ONR  Task  No.  NR  334-003. 
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V • n - u • n 


3#-  3<SfU3^s  , 3^fn^s  - t , (1.4) 

V • n is  the  outward  normal  velocity  of  the  flow  region  6^(t),  and  Vg  . n is 
prescribed  on  30,(0 . 3^  is  called  the  "free  boundary".  (For  a more  detailed 

discussion  of  the  asymptotic  conditions  (1.3a),  see  chapter  II  of  the  report  by 
Rogers  (Ref.  6).) 

For  classical  solutions  of  the  differential  equations  (1.1),  the  velocity 
field  u(x,t)  will  possess  derivatives  with  respect  to  x and  t.  When  an  initially 
differentiable  velocity  field  Uq  evolves  into  a field  u(x,t)  which  is  no  longer 
differentiable,  it  becomes  necessary  to  re-interpret  the  problem  (1.1)  - (1.4) 
in  a suitably  generalised  manner.  Such,  in  fact,  is  the  case  when  a wave  spills 
over  and  falls  back  on  the  surface.  At  the  moment  of  impact  with  the  surface, 
the  velocity  field  at  the  point  of  Impact  is  discontinuous. 

Pursuing  further  the  problem  of  a wave  spilling  over,  or  in  general  of 
the  collision  of  two  incompressible  fluids  with  free  surfaces,  we  may  consider 
the  idealized  problem  shown  in  Figure  1 . Here  we  imagine  that  we  have  a cylind- 
er in  a region  of  space  where  there  is  no  gravity,  and  that  this  cylinder  con- 
tains two  equal  masses  of  liquid  which  are  moving  without  friction  with  equal 
and  opposite  velocities  along  the  axis.  The  free  surface  of  each  liquid  mass 

consists  of  two  components  orthogonal  to  the  cylinder  axis.  At  the  moment  of 

% 

collision  we  expect  the  condition  (1.1b)  on  the  velocity  to  be  violated.  In- 
stead, if  we  denote  distance  along  the  cylinder  axis  by  z,  set  the  origin  in  the 
plane  of  collision,  and  denote  the  speed  of  each  liquid  mass  before  collision  by 
U,  we  get  at  the  moment  of  collision 

7 • u - - 2 U6(z)  . (1.5) 

Thus  this  example  indicates  that  we  will  need  to  re-evaluate  the  condition  (1.1b) 
in  expressing  the  laws  of  hydrodynamics  in  a form  suitable  for  treatment  of  the 
general  free  boundary  problem. 

A different  sort  of  problem  arises  in  connection  with  the  delineation  of 

the  free  boundary.  When  equations  (1.1)  - (1.4)  have  a classical  solution,  one 

finds  from  the  requirement  of  consistency  between  (1.1a)  and  (1.1b)  that 

3u.  3u, 

AP " * po  £ ^ ^7  (1-6> 

The  need  for  consistency  between  (1.1a)  and  (1.3b)  leads  to  a Neumann  condition 
on  the  pressure  at  3(^# . This,  the  asymptotic  condition  (1.3a),  and  the  free 
surface  condition  (1.3c)  combine  with  (1.6)  to  determine  the  pressure  throughout 
Cl.  With  P determined,  (1.1a)  can  be  solved  for  u,  and  (1.3d)  then  serves  to 
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determine  the  time  evolution  of  </^(t).  In  the  case  when  • 0,  we  find  that 
the  problem  (1.1)  - (1.4)  is  invariant  under  the  transformation 

x'  - cot,  t'  - a*t  (1.7) 

for  any  constant  a > 0.  What  this  means  is,  for  example,  that  if  a flow  start- 
ing from  rest  with  a particular  initial  surface  deformation  results  in  the  fall- 
ing over  and  splashing  down  of  a wave  in  time  Tq,  then  the  problem  with  similar 
Initial  conditions,  except  that  the  surface  deformation  is  scaled  relative  to 
that  in  the  first  problem  by  a factor  a < 1 , results  in  the  falling  over  and 
splashing  down  of  a wave  in  time  a*4^.  Accordingly,  we  may  generally  expect 
that,  for  any  flow,  small  perturbations  on  the  flow  will  have  the  property  that 
waves  on  sufficiently  small  scales  will  be  continually  falling  over,  and  that 
the  topology  of  such  perturbed  flows  may  not  even  be  determinate.  In  the  case 
of  the  breaking  of  a single  wave,  but  even  more  so  in  this  case  of  multiple 
breaking,  the  determination  of  the  evolution  of  the  free  surface  through  an 
equation  like  (1.3d)  becomes  ambiguous.  For,  if  we  regard  (1.3d)  as  an  ordinary 
differential  equation  for  the  motion  of  a fluid  element  on  the  free  surface, 
such  an  equation  does  not  have  a unique  solution  when  the  velocity  field  u is 
discontinuous  at  the  free  surface.  In  addition,  it  is  open  to  question  whether 
the  portion  of  the  free  boundary  contained  in  any  unit  ball  will  have  bounded 
measure,  or  whether  the  points  on  the  free  boundary  will  be  regular  points  for 
the  Poisson  equation  (1.6).  If  they  are  not  regular  points,  the  meaning  of  the 
condition  (1.3c)  will  need  re-examlnatlon. 

2.  Role  of  Conservation  Laws 

Experience  with  other  free  boundary  problems  has  shown  that  the  equation 
which  drives  the  motion  of  the  free  surface  is  often  merely  an  expression  of  a 
fundamental  conservation  law.  For  example,  in  the  Stefan  problem,  conservation 
of  energy  is  paramount  (Ref.  1)  and  in  some  model  hyperbolic  problems,  the  shock 
conditions  are  an  expression  of  other  conservation  laws  (Ref.  7).  Recognition 
of  this  fact  has  made  it  possible  to  devise  topology-independent  algorithms  to 
solve  such  free  boundary  problems.  The  use  of  such  methods  is  especially  de- 
sirable in  the  present  problem  where,  as  we  have  indicated,  the  surface  may  be- 
come so  complicated  that  it  may  become  impossible  to  follow  its  motion,  not  only 
in  practice,  but  also  in  theory. 

The  question  of  what  conservation  laws  are  appropriate  in  order  that  our 
mathematical  model  adequately  represent  a physical  situation  is,  of  course,  a 
problem  of  physics.  In  our  relatively  simple  situation,  the  answer  seems  rather 
clear.  The  best-known  conservation  laws  of  Newtonian  physics  are  those  of  mass, 
momentum,  and  energy.  For  the  classical  flow  of  an  lnvlscld  liquid  one  may  de- 
rive energy  conservation  from  the  conservation  laws  for  mass  and  momentum  and 
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generally  energy  conservation  plays  a subsidiary  role,  being  derivable  from  the 
equations  of  motion  in  reversible  physical  situations  and  requiring  reformula- 
tion in  Irreversible  situations.  In  hydrodynamic  theory  the  precedence  of  mass 
and  momentum  conservation  over  energy  conservation  is  assumed,  for  example,  in 
the  derivation  of  the  jump  conditions  for  solutions  of  the  nonlinear  shallow 
water  equations  (Ref.  12).  The  energy  which  is  lost  is  assumed  somehow  dissi- 
pated in  other,  irreversible  processes,  or  in  turbulence.  This  relation  of 
energy  loss  to  irreversibility  is  a natural  complement  to  the  connection  of 
energy  conservation  with  temporal  homogeneity  in  Hamiltonian  mechanics.  (In 
dynamical  systems  which  are  richer  in  degrees  of  freedom  than  ours,  the  burden 
of  irreversibility  is  shifted  from  the  energy  to  the  entropy.) 

Regarding  the  conservation  laws,  one  notes  that  (1.1a)  is  a statement  of 
conservation  of  momentum  (in  the  case  g - 0) . For  a fluid  whose  elements  do  not 
undergo  a density  change  as  they  move  and  whose  velocity  is  uniformly  differen- 
tiable from  one  point  to  another,  so  that  the  trajectories  of  different  elements 
remain  distinct,  (1.1b)  is  a statement  of  mass  conservation.  When  the  density 
varies  discontinuously  in  space,  as  it  does  at  the  water  surface,  the  governing 
equation  for  the  free  surface,  (1.3d),  is  also  an  expression  of  mass  conserva- 
tion. Thus,  it  would  appear  that  at  least  one  of  the  problems  referred  to  above, 
how  the  free  surface  evolves  in  time,  may  be  resolved  by  writing  a law  of  mass 
conservation.  When  the  density  p and  velocity  u are  differentiable,  this  takes 
the  form 

Pt  + V • (pu)  - 0 (2.1) 

As  we  have  noted,  classically  one  may  think  of  (1.1b)  as  an  "equation  of 
state"  in  terms  of  which  the  pressure  is  determined.  In  the  classical  picture 
the  equation  of  state  is  a constraint,  and  in  the  process  of  satisfying  this 
constraint  the  momentum  is  altered  by  the  term  -VP.  We  go  one  step  further,  and 
suggest  that  in  the  absence  of  a constraint  the  pressure  vanishes,  that  is,  the 
pressure  arises  only  when  the  constraint  cannot  be  satisfied  without  it.  In 
accordance  with  the  kinetic  theory  interpretation  of  pressure,  we  may  identify 
the  pressure  with  the  transfer  of  momentum  across  the  surface  of  a fluid  element 
in  the  direction  of  its  normal  brought  about  by  the  action  of  the  constraint. 

The  boundary  condition  (1.3c)  suggests  that  at  the  free  boundary  the  constraint 
is  automatically  satisfied. 

Let  us  inquire  further  what  the  nature  of  the  constraint  should  be  for 
the  generalized  flows  of  interest  to  us.  To  this  end,  we  reconsider  the  situa- 
tion depicted  in  Figure  1.  We  make  two  observations:  First,  the  dynamics  of 
the  system  should  in  no  way  be  affected  by  the  way  we  extend  the  velocity  field 
to  the  region  where  p ■ 0 — it  is  momentum  and  not  velocity  which  is  dynamically 
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important.  Second,  there  Is  no  problem  with  a velocity  distribution  whose  di- 
vergence approaches  (l.S)  at  the  moment  of  collision.  The  fundamental  fact  re- 
garding the  collision  that  we  want  to  maintain  is  that  the  two  bodies  of  fluid 
should  not  enter  each  other.  That  is,  the  fluid  density  should  not  exceed  pg. 

We  write  this  as  a one-sided  constraint 

P < P0-  (2.2) 

We  will  see  that  (2.2)  is  the  appropriate  generalization  of  the  classical  con- 
straint (1.1b)  and  also  of  the  boundary  condition  (1.3c).  Its  one-sided  nature 
reminds  us  of  variational  inequalities. 

3.  Algorithmic  Description  of  Hydrodynamics 

One  of  our  governing  equations  is  the  conservation  of  mass,  (2.1).  In 
accordance  with  the  observations  of  the  last  section,  in  the  absence  of  the  con- 
straint we  will  have  the  equations  of  momentum  "conservation": 

(pu)t  + V • (puu)  - - pg  £ (3.1) 

Equations  (3.1)  and  (2.1)  form  a set  of  hyperbolic  conservation  laws.  For  the 
actual  hydrodynamic  flow,  we  will  solve  them  subject  to  the  constraint  (2.2). 

Of  course  (2.1),  (3.1),  and  (2.2)  are  generally  Inconsistent,  and  we  have 
to  make  clear  what  we  mean  by  "solving"  (2.1)  and  (3.1)  subject  to  (2.2).  We 
will  do  this  by  giving  an  algorithm  for  approximating  the  solution  of  the  evo- 
lutionary problem  in  which  we  start  with  initial  data  p(0),  u(0)  and  try  to  find 
p(t) , u(t)  for  t > 0.  Our  algorithm  will  be  dependent  on  a parameter  T which 
we  call  the  time  step,  and  will  generate  from  a pair  of  quantities  p,  pu,  with  p 
satisfying  (2.2),  another  pair  of  quantities  p,  pu,  with  p satisfying  (2.2). 

We  denote  the  result  of  this  operation  symbolically  as 


(p,  pu)  - S(t)  (p,  Pu)  . (3.2) 

When  we  speak  of  "solution"  of  the  problem,  we  mean  that  for  t > 0 the  operators 


(s  (£))”  - S*(t) 


(3.3) 


as  n -*■  <*>.  (3.3)  is  to  be  understood  to  hold  in  an  appropriate  function  space. 

More  will  be  said  about  this  in  our  next  lecture  (Ref.  9).  However,  it  be- 
hooves us  to  point  out  that  we  have  not  yet  proven  the  crucial  step  (3.3),  and 
thus  we  cannot  speak  of  a "solution"  of  the  problem  in  any  rigorous  mathematical 
sense.  In  this  lecture  we  content  ourselves  with  an  indication  that,  when  the 
flow  quantities  have  sufficient  regularity  in  space  and  time,  our  algorithm  re- 
duces to  an  approximate  algorithm  for  solving  the  Euler  equations,  which  may  be 
expected  to  converge  to  the  actual  solution  as  T 0,  under  the  same  presupposi- 
tions regarding  regularity.  The  next  lecture  will  focus  on  the  sorts  of 
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solutions  we  expect  to  emerge  from  the  analysis  of  convergence,  and  the  sense  In 
which  the  lnvlscid  hydrodynamic  initial  value  problem  may  generally  be  regarded 
as  well-posed.  But  the  problems  of  convergence  and  regularity  of  the  flows  con- 
verged to  are  ongoing  problems,  presently  uncompleted. 

Our  algorithm  "solves"  (2.1)  and  (3.1)  subject  to  (2.2)  in  the  following 
sense:  The  hyperbolic  conservation  laws  (2.1)  and  (3.1)  are  "solved"  for  a time 
interval  T.  Then  the  densities  of  mass  and  momentum  are  adjusted  to  satisfy 
(2.2),  in  a manner  consistent  with  global  conservation  of  mass  and  momentum. 

When  we  say  that  incompressible  flows  may  be  considered  to  evolve  through  a 
system  of  conservation  laws  in  conjunction  with  a constraint,  our  statement  is 
premised  on  the  conjectured,  but  as  yet  unproven,  existence  of  the  limit  of  the 
algorithm  as  T -*■  0.  There  is  no  doubt  a certain  lack  of  elegance  in  our 
approach  through  a family  of  solutions  dependent  on  a parameter  T,  but  it  is 
perhaps  no  worse  than  the  situation  which  arises  in  making  precise  the  solution 
of  an  initial  value  problem  for  an  ordinary  differential  equation. 

In  what  follows,  we  will  attempt  a reasonably  complete  description  of  the 
algorithm  which  is  to  render  an  approximation  to  the  flow.  However,  what  we 
present  is  by  no  means  our  first  approach  to  the  problem,  and  along  the  way 
mathematical  simplifications  have  arisen  which  have  removed  the  algorithm  some- 
what from  its  pristine  physical  orientation.  For  a more  complete  description  of 
the  physical  considerations  which  led  us  to  make  some  of  our  initial  choices,  we 
refer  the  reader  to  a more  complete  write-up  (Ref.  10).  What  we  present  here  is 
a mathematical  object,  which  will  rise  to  the  status  of  theory  or  fall  into 
disrepute  according  to  its  internal  consistency.  No  doubt  later  versions  will 
differ  in  detail,  but  we  suspect  that  the  main  elements  will  remain  Intact. 

To  "solve"  the  conservation  laws  (2.1)  and  (3.1)  for  a time  r,  we  intro- 
duce a distribution  function  F(x,v,t)  satisfying  the  collisionless  Boltzmann 
equation 

Ft+v-VF-g|^-0  , 0 < t < T , (3.4a) 

z 

and  initial  conditions 


F(x,v,0)  - p(x)  6 (v-u(x))  . (3.4b) 

It  is  easiest  to  give  boundary  conditions  for  F in  terms  of  the  characteristics, 
whose  equations,  away  from  boundaries,  are 


dx  dv 

dt  " v>  dt 


gk 


(3.5) 


These  are  just  the  equations  of  classical  particles  moving  under  the  Influence 
of  gravity  without  collisions.  At  a rigid  boundary,  we  require  that  the  equa- 
tions of  the  characteristics  describe  the  trajectories  of  particles  reflecting 
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specularly  from  the  boundary.  Thus,  If  a characteristic  strikes  a rigid  boundary 
at  time  t£  and  the  velocity  of  the  rigid  boundary  is  Vg  at  the  point  and  time 
where  It  Is  struck  by  the  characteristic,  we  set 

x(t*)  - x(t")  , (3.6a) 

c c 

v(t*)  - nx  (v(t‘)xn)  + (2V  • n - v(t“)  • n)n,  (3.6b) 

c c sc 

where  n is  the  unit  outward  normal  to  the  boundary  at  the  point  and  time  referr- 
ed to. 

Finally,  we  determine  approximate  solutions  of  (2.1)  and  (3.1)  at  time  T 
through 

p - /F(x,v,t)dv  f (3.7a) 

f\y\j 

pu  - /vF(x,vtT)dv  . (3.7b) 

In  physical  language,  equations  (3.7)  state  that  all  fluid  elements  at  the 
same  location  after  the  passage  of  time  T have  collided  inelastlcally.  We  point 
out  that  this  assumption  that  collisions  are  inelastic  is  not  mandatory,  but  it 
seems  like  a simple  and  reasonable  first  approximation  for  the  problems  that  in- 
terest us.  Other  assumptions  are  possible.  The  allowance  of  inelastic  colli- 
sions permits  the  decay  of  energy,  and  an  element  of  irreversibility  enters  into 
our  algorithm,  although  the  Euler  equations  themselves  are  formally  reversible 
in  time.  We  note  that  some  additional  assumptions  regarding  the  nature  of  colli- 
sions have  been  needed  to  make  the  evolutionary  problem  determinate  in  the  gen- 
eral case,  and  our  treatment  of  the  conservation  laws  has  provided  a set  of  such 
assumptions.  For  example,  in  the  situation  depicted  in  Figure  1,  a number  of 
possibilities  after  collision  will  be  consistent  with  the  requirements  we  have 
made  heretofore.  One  possibility  is  for  the  two  liquid  masses  to  collide  and 
them  come  to  rest  instantaneously,  with  all  energy  lost  inelastlcally  at  the 
moment  of  impact.  Another  possibility  is  for  them  to  collide  totally  elastically, 
bouncing  off  one  another,  with  the  flow  totally  reversible.  There  are  also  inter- 
mediate possibilities,  with  a loss  of  speed  for  all  the  fluid  being  one,  and 
with  some  of  the  fluid  being  brought  to  rest  and  the  remainder  rebounding  elas- 
tically being  another.  Although  we  are  getting  somewhat  out  of  sequence,  since 
we  have  not  described  how  the  algorithm  treats  the  constraint  condition  (2.2) 
yet,  we  note  that,  according  to  the  assumption  of  inelastic  collisions  made  in 
(3.7),  in  the  limit  as  T 4 0 for  the  case  shown  in  Figure  1,  we  will  get  the 
first  possibility  listed  above. 

Classical  flows  in  which  the  velocity  is  Lipschitz  continuous  in  space  un- 
iformly in  time  will  not  permit  the  collision  of  fluid  elements  for  T sufficiently 
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small,  and  thus  this  element  of  irreversibility  will  not  enter.  We  will  refer 
to  the  conservation  laws  (2.1)  and  (3.1)  with  g - 0 as  the  higher-dimensional 
form  of  Burgers'  equation.  We  have  noted  elsewhere  (Ref.  7)  how  in  one  dimension 
the  proper  solution  of  (2.1)  and  (3.1),  as  outlined  in  (3.4)  and  (3.7),  differs 
from  the  solution  of  the  formally  equivalent  conservation  law 


We  will  say  nothing  in  this  lecture  about  the  convergence  as  T ♦ 0 of  the  solu- 
tion of  the  conservation  laws  outlined  in  (3.4)  and  (3.7).  The  subject  will  be 
raised  in  the  next  lecture. 

Assuming  that  we  have  adequately  solved  (2.1)  and  (3.1),  let  us  see  how 

<\, 

we  satisfy  the  constraint  (2.2).  If  p _<  p , the  constraint  has  no  effect,  and 

_ <\j  — vv*  ^ ” 

we  set  p ■ p,  pu  ■ pu.  If  p > pQ  somewhere,  in  violation  of  (2.2),  we  have  to 
realize  that  during  the  time  T while  mass  was  being  convected  according  to  (3.4), 
(3.5),  (3.6),  and  (3.7a),  other  processes  were  also  taking  place.  It  may  help 

f\j 

to  observe  that  p is  a linear  functional  of  p,  and  that  we  may  envisage  p as 
the  accumulation  of  Independently  moving  mass  densities,  or  "streams".  The 
other  processes  that  took  place  in  the  time  T were  of  the  following  sort: 

Whenever  there  was  an  accumulation  of  mass  yielding  a density  > Pq,  the  particles 
in  the  region  of  excess  density  were  considered  to  be  undergoing  rapid  elastic 
collisions  which  resulted  in  their  spreading  out  from  the  region  of  density  ex- 
cess in  an  isotropic  manner.  As  new  density  excesses  arrived  at  a point  from 
additional  streams,  they  also,  in  addit^.  to  the  excess  particles  which  had  not 
yet  spread  out  from  previous  collisions  and  which  therefore  still  contributed 
to  the  density  excess  at  the  point,  underwent  such  elastic  collisions  with  a 
resultant  spreading  out.  Such  collisions  occurred  with  extreme  rapidity,  with 
the  result  that  after  the  time  t all  streams  which  had  contributed  to  the  den- 
sity excess  at  a point  had  spread  out  and  no  excess  was  left. 

Each  group  of  collisions  with  attendant  spreading  out  of  mass  was  iso- 
tropic, and  there  were  many  such  processes  going  on,  until  a sort  of  "steady 
state"  was  achieved.  Now  in  fact  the  effect  of. an  isotropic  spreading  out  of  a 
mass  distribution  do(x)  is  to  replace  do  by 

(/P(x,x’)do(x'))dx 

where 

/P(x,x')dx  - I , (3.8a) 


/P(x,x') (x-x')dx  - 0 


(3.8b) 
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/P(x,x') (x-x'Jj^x-x') jdx  » q(x')61j  (3.8c) 

Repeated  application  of  minute  (q  •*  0)  processes  of  this  sort  is  equivalent  to 
repeated  application  of  Gaussian  distributions  (Ref.  4),  and  thus  we  may  replace 
P by  a Gaussian  with  small  variance: 


P(x,x')  -*■ 


(x-x')^/4da(x') 


(4irda(x' 


jn 


(3.9) 


where  N is  the  dimensionality  of  the  space.  (3.9)  holds  in  the  interior  of  the 
fluid.  Of  course,  at  a rigid  boundary  any  fluid  which  spreads  out  cannot  pass 
through  the  boundary.  Instead,  there  must  be  a reflection  at  the  boundary,  as 
there  was  in  the  case  described  by  equations  (3.6).  So  in  general,  we  replace 
(3.9)  by 


P(x,x')  -*■  e^  da^x  ^ - S+(da(x' 


)) 


(3.10) 


where  S (da(x'))  is  the  semigroup  generated  by  the  Laplace  operator  for  the  re- 
gion exclusive  of  rigid  bodies,  with  the  requirement  of  zero  normal  gradient  at 
the  rigid  boundaries. 

In  the  hydrodynamic  case  only  that  part  of  the  mass  distribution  corre- 
sponding to  density  exceeding  pQ  spreads  out,  and  thus  the  operator  S+(da(x')) 
acts  only  on 


where 


f(p)  - 


f(p(x'))dx’ 


lp-po  p±p0 

I 

‘ 0 p 1 po 


(3.11) 


and  p is  a general  mass  density.  What  happens,  then,  is  that  an  Initial  mass 

A 

density  p is  replaced  by 


P(1)  - p - f(p)  + S+(da1(x,)f(p)  = Fjp  , 


(3.12a) 


^ ( 1 ) 

p is  replaced  by 


P(2)  - F2£(1)  - p(1)  - f(p(1))+  S+(da2(x’))  f(p(1))  , (3.12b) 

A (fO 

and  in  general  p is  replaced  by 

P(n+1)  - Fn+1  p(n)-  p(n)  -f(p(n))  + S+(dan+1(x'))  f (p(n)).  (3.12c) 

After  many  such  collisions  n + »»e  achieve  a steady  state. 


A’ 


<**■ 


- 10  - 


J 

j 


Referring  to  (3.10)  end  letting  da(x')  0,  we  see  that  the  steady  state 
Is  the  steady  state  (a  «*>)  of  the  equation 

6a  - A+  (y(x,a)  f (§))  , (3.13a) 

0 (a  - 0)  - p (3.13b) 

where 

U(x,a)  - , da  - 8“p  da(x,a)  , (3.13c) 

and  the  assumption  that  collisions  take  place  wherever  there  is  a mass  density 
excess  is  reflected  In  the  condition 


y(x,a)  > 0 (3.13d) 

p has  been  envisioned  as  the  accumulation  of  a number  of  Independently  moving 
streams,  and  hence  to  find  the  new  mass  density  after  all  the  mass  redistribu- 
tions due  to  collisions  have  taken  place,  we  should  solve  equation  (3.13)  with 

/s 

an  initial  density  p and  an  inhomogeneous  term  on  the  right-hand  side  reflecting 
the  addition  of  other  contributions  to  as  the  parameter  a runs  from  0 to  00  . 
However,  we  have  found  (Ref.  8)  that  the  steady  state  is  independent  of  the 
order  in  which  contributions  to  p are  inserted  into  the  equation.  A reduction 
in  y(x,a)  may  be  viewed  as  a change  in  the  order  in  which  contributions  appear, 
and  thus  we  note  that  the  steady  state  is  the  same  as  that  for  the  problem 

0a-A+f(e),  (3.14a) 

0(a-  0)  - p . (3.14b) 

(3.14)  is  recognized  as  a one-phase  Stefan  problem.  In  terms  of  the  solution 
of  this  problem,  the  new  mass  density  is 

5 - 6 (<»)•  (3-15) 

(3.14)  and  (3.15)  serve  to  determine  the  location  of  the  hydrodynamic  free 
boundary. 

We  turn  next  to  the  effect  of  the  mass  redistribution  on  the  momentum  den- 
sity. Just  as  the  elastically  colliding  particles  carry  a mass  with  them  as 
they  move,  they  also  carry  a velocity  u,  which  is  as  yet  undetermined.  But  in 
addition,  since  velocity  is  nothing  but  rate  of  spatial  displacement,  the 
particles  must  have  associated  with  them  a momentum  due  to  the  fact  of  their  re- 
distribution. Since  all  these  processes  take  place  in  a time  t,  to  lowest 
order  in  x we  may  associate  with  a particle  which  has  moved  from  x’  to  x the 
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it: 
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velocity  (x-x')/t.  Away  from  a rigid  boundary,  in  the  process  represented  by 
(3.12a),  the  momentum  density  pu  will  be  replaced  by 

(PU)V  ' - Pu  - U f(p)  + S(dal(x,))(f(P)u) 


i 


(4irda1(x')J,/2  T 


x-x'  e“(x-x,)  /4dal(x,)f(S(x'))dx’ 


AA  ^ 

- pu  + (S^ot^x'))-!)  (G  f(p)) 

V (S(dai(x'))  f^x’^da^x'))  , (3.16a) 

where 

S(da1(x'))  - eAd0tl(x,)  (3.16b) 

The  process  represented  by  (3.12c)  will  result  in  the  replacement  of  (pu) ^ by 
(pG)(n+1)  - (pG)(n)  + (S(dan+1(x'))  - 1)  (u  f(P(n))) 

- \ V(S(dan+1(x’))f(G(n)(x'))lan+1(x'))  . (3.16c) 


Letting  n -*•  <»  and  da^x')-*  0,  we  get,  independent  of  the  order  in  which  the 
colllsional  processes  associated  with  the  mass  redistribution  occur, 

pu  - pu  + A(uv)  - ^ Vv  (3.17a) 

where 

00 

v “Jo  f(0>da  (3.17b) 

and  6 satisfies  (3.14).  (3.17)  is  to  be  solved  subject  to  the  boundary  condi- 

tion (3.6)  on  u at  the  rigid  boundary  3^fi(t).  (More  precisely,  the  normal  com- 
ponent u • n satisfies  (1.3b),  and  the  derivative  in  the  normal  direction  of 
the  tangential  component  uXn  vanishes.) 

It  may  seem  that  there  is  some  mystery  associated  with  our  determination 
of  a velocity  field  in  terms  of  the  displacement  of  moving  particles  in  a time 
interval,  as  opposed  to  its  determination  through  higher  order  time  derivatives 
of  the  displacement,  namely,  the  acceleration.  However,  the  determination  here 
really  has  grown  out  of  a consistency  argument,  and  we  can  give  an  example  from 
elementary  mechanics  to  illustrate  our  point. 

Consider  the  situation  shown  in  Figure  2.  We  have  a particle  moving  on 
the  surface  of  a rigid  body  under  the  Influence  of  gravity.  The  particle  may 
move  on  the  surface  or  above  it,  but  may  not  enter  the  rigid  body.  Thus  we 
have  a one-sided  constraint  on  the  motion  of  the  particle,  similar  in  some 
respects  to  (2.2).  One  may  devise  an  algorithm  to  determine  the  motion  of  the 
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particle  as  follows:  Given  the  particle  position  and  velocity  (r,u)  at  a given 
time,  we  let  the  particle  follow  the  familiar  parabolic  path  appropriate  to  mo- 
tion in  a gravitational  field  without  any  constraints.  This  carries  the 
particle  to  P after  a time  T,  when  it  has  velocity  u.  If  P lies  on  or  above  the 
rigid  body,  no  correction  is  necessary,  and  we  can  set  P ■ If,  u ■ u for  the  new 
position  and  velocity  of  the  constrained  motion.  If  P lies  inside  the  rigid 
body,  we  satisfy  the  constraint  by  moving  the  particle  back  to  the  nearest  point 
P on  the  body's  surface.  Then,  to  be  consistent  with  the  fact  that  velocity  is 
rate  of  spatial  displacement,  we  have  to  add  to  u the  displacement  from  ^ to  P 

r\j 

divided  by  T.  This  addition  to  u is  known  in  mechanics  as  the  normal  force 
exerted  on  the  particle  by  the  body.  At  a point  Pf  on  the  body  the  particle  may 
leave  the  surface.  Pf  may  be  thought  of  as  a "free  boundary". 

We  regard  the  region  of  flow  where  0 < p < to  be  a "spray".  This  is 

more  a mathematical  artifice  than  a physically  complete  representation  of  an 
actual  spray.  A more  detailed  description  of  some  of  the  physical  assumptions 
made  in  our  characterization  of  the  fluid  in  the  region  where  0 < p < as 
a spray  is  given  elsewhere  (Ref . 10) . The  possibility  of  the  development  of  a 
spray  in  the  non-classical  formulation  of  hydrodynamics  is  analogous  to  the 
possibility  of  "slush"  formation  in  the  non-classical  formulation  of  the  Stefan 
problem  (Ref.  I).  Indeed,  as  seen  in  (3.11)  and  (3.14),  there  is  a clear  cor- 
respondence between  the  enthalpy  and  latent  heat  in  the  one-phase  Stefan  pro- 
blem, on  the  one  hand,  and  the  mass  and  liquid  density  in  hydrodynamics,  on  the 
other.  Similar  interpretations  may  likewise  be  given  to  "npray"  and  "slush". 

In  the  latter  case,  slush  occupying  a region  of  positive  measure  may  be  con- 
ceived as  a mixture  of  minute  volumes  of  two  phases  of  a substance,  such  that 
the  volume  of  each  phase  has  a positive  measure  in  each  subset  of  positive 
measure  in(^  . In  the  former  case,  we  think  of  spray  occupying  a region  (?t  of 
positive  measure  as  consisting  of  minute  volumes  of  liquid  (p  • Pq)  and  vacuum 
(p  - 0) , with  the  volume  of  each  phase  in  each  subset  of  0^  of  positive  measure 
having  positive  measure.  As  yet  we  do  not  have  any  examples  of  flows  in  which 
we  can  show  rigorously  that  sprays  must  develop  in  order  for  a solution  of  the 
equations  to  exist.  Nevertheless,  as  we  pointed  out  in  the  Introduction,  it  is 
by  no  means  clear  that  the  hydrodynamic  free  boundary  can  always  be  sharply  de- 
fined, and  we  leave  open  the  possibility  of  the  formation  of  a diffuse  free 
boundary. 

In  concluding  this  section  dealing  with  an  algorithmic  representation  of 
a generalized  hydrodynamics,  we  remark  that  numerical  results  based  on  the  al- 
gorithm have  been  obtained,  and  are  currently  being  prepared  for  publication 
(Ref.  11).  The  numerical  treatment  of  the  hyperbolic  conservation  laws  follows 
the  path  laid  down  in  (3.4)  - (3.7).  The  steady  state  of  the  one-phase  Stefan 


1-4 ; , __ 
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problem  (3.14)  Is  found  using  the  algorithm  (3.12)  for  the  case  da^x')  - da 
(Ref.  1).  The  linear  elliptic  problem  (3.17)  Is  solved  by  finding  the  steady 
state  of  a parabolic  equation,  which  Is  in  turn  solved  through  a variation  on 
an  algorithm  applicable  to  a class  of  hyperbolic  and  parabolic  problems  (Ref.  2). 
4.  Consistency  of  the  Algorithm 

To  demonstrate  consistency  of  the  reformulated  hydrodynamics  with  clas- 
sical hydrodynamics  in  the  regime  where  the  latter  is  meaningful,  we  examine 
alternate  forms  of  our  equations  when  the  velocity  field  is  differentiable. 
Consistency  will  be  demonstrated  if  we  show  that  the  lowest  order  terms  in  T 
are  identical  In  both  formulations.  The  consistency  of  (3.4)  - (3.7)  with  the 
hyperbolic  conservation  laws  (2.1)  and  (3.1)  is  quite  straightforward.  As  we 
have  pointed  out,  when  the  velocities  are  differentiable  no  collisions  of  fluid 
elements  will  occur  for  T sufficiently  small,  and  (3.4)  - (3.7)  will  then  solve 
(2.1)  and  (3.1)  exactly  in  the  Interior  of  the  flow  region.  As  regards  bound- 
ary conditions,  it  follows  from  (3.6)  that  the  average,  over  a small  time  in- 
terval, of  the  component  of  momentum  normal  to  a rigid  boundary  must  approach 
that  component  of  velocity  of  the  boundary  times  the  average  over  the  same  time 
interval  of  the  density,  as  the  distance  to  the  boundary  approaches  zero.  This 
is  consistent  with  (1.3b). 

Let  us  then  focus  attention  on  the  second  half  of  the  algorithm,  which 
deals  with  the  ramifications  of  the  constraint  (2.2).  If  we  can  show  that  (3.17) 
goes  over,  to  first  order  in  T,  to  the  equation 

pu  - pu  - TVP  (4.1) 

where  P satisfies  (1.6),  we  will  have  shown  that  our  algorithm  reverts  to  a 
split-step  scheme  for  solving  the  classical  (1.1).  Comparing  (4.1)  with  (3.17), 
we  see  that  the  obvious  correspondence  to  make  is  that 

^♦P  (4.2) 

T 

(3.17b),  (3.14),  and  (3.15)  lead  to 

A+  v - p - p (4.3) 

Suppose  at  a given  time  we  have  a density  P f.  Pq  and  a velocity  u satisfying 
(1.1b).  In  the  interior  of  the  liquid,  p«  Pq.  Now,  assume  the  velocity  is 
differentiable,  so  that  (2.1)  and  (3.1)  imply 

+ u . Vu  ■ -gk  (4.4) 

% 

Integrating  (4.4)  over  the  time  interval  T to  get  u,  we  find  to  flret  order 

I * 
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in  T 


i '■t' 

» » ‘ 

Vk/  - 


3u,  3u. 

The  equation  of  mass  conservation. 


(4.5) 


Pt  + u 


Vp  - - p7 


gives  for  p in  the  interior  of  the  region  p*pn,  to  second  order  in  T 


P “ Pr 


0* 

P0  I * u + V • “) 


t2  r5»i  V) 
■p0  + p0^^  < 


(4.6) 


For  points  Interior  to  the  liquid  region  at  the  given  time  (p  ■ Pq)  and  also 
the  liquid  region  a time  step  later  (p  • Pq) , it  follows  from  (4.3)  that 


p0  T‘  ~^i 

2 **  3xj  3x1  * 


Av  ■ - 


in  agreement  with  (1.6)  and  the  correspondence  (4.2).  We  recall  that  (1.6)  is 
Just  the  condition  to  make  V • u ■ 0 to  lowest  order  in  x. 

With  respect  to  boundary  conditions,  if  we  have  a classical  flow  with  a 

f\j 

sharp  free  boundary,  p will  fall  rapidly  from  the  expression  (4.6)  to  0 at  the 
free  boundary.  In  the  Interior  of  the  liquid,  as  x •+■  0,  we  find  from  (4.6) 
and  (4.3)  that  |Av|  is  small  compared  to  Pq.  whereas  outside  the  free  boundary 
we  will  have  approximately  Av  - Pq.  Thus,  although  one  derives  from  (3.17b) 
and  (4.3)  that  v » Vv  ■ 0 at  the  free  boundary,  there  will  be  a sort  of  "bound- 
ary layer"  there  in  which  Vv  changes  from  0 to  a finite  value.  Since  by  (4.6) 

«\,  2 

p will  differ  from  P0  by  0 (x  ) in  the  interior  of  the  liquid,  mass  conserve- 

^ 2 
tion  will  require  that  this  "boundary  layer"  have  thickness  0(x  ) if  the  free 

2 

boundary  has  bounded  curvatures.  From  (4.3),  Vv  will  change  by  0(x  ) over  this 

4 

boundary  layer  and  v will  change  by  0(x  ).  Just  inside  the  boundary  layer 
/2/x^]Vv  will  assume  a value  which  does  not  necessarily  vanish  as  x 0,  but 
|2/x^J  v will  «*•  0 as  T 0.  Hence  to  lowest  order  in  X the  boundary  condition 
(1.3c)  will  be  redeemed. 

Something  similar  occurs  at  the  rigid  boundaries.  On  account  of  (3.14), 
we  get  n • Vv  » 0 at  rigid  boundaries.  However,  (3.4)-(3.7)  predict  that, 

over  a "boundary  layer"  with  thickness  0(gx2) , p - pn  will  be  0(p  ).  Across 

2 u 0 

this  boundary  layer  2/x  n • Vv  will  Jump  from  0 to  a value  O(pQg).  Interpret- 
ing the  asymptotic  condition  (1.3a)  as  an  approximation  to  the  case  where  the 
fluid  is  bounded  below  by  a portion  of  a rigid  plane  situated  at  a large  nega- 
tive value  of  z,  we  easily  confirm  the  validity  with  respect  to  this  condition 
of  the  correspondence  (4.2).  Similar  results  obtain  at  other  rigid  boundaries, 
but  we  note  that  agreement  between  the  different  formulations  is  built  in  by 
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requiring  solutions  of  (3.17)  to  have  the  derivative  in  the  normal  direction  of 
their  tangential  component  vanish,  and  to  have  the  normal  component  satisfy 
(1.3b). 

With  the  correspondence  (4.2),  the  term  A(uv)  on  the  right  of  (3.17a)  is 
2 

seen  to  be  0(t  ),  and  thus  to  have  no  effect  on  the  consistency  of  our  formula- 
tion with  the  classical  equations  (1.1)  - (1.4).  Reference  to  (3.16)  and  the 
discussion  preceding  (3.16)  shows  that  the  extra  term  represents  the  fact  that 
the  elastically  colliding  particles  in  our  picture  carry  the  mean  velocity  u 
with  them.  Even  though  this  term  may  be  removed  from  (3.17a)  without  changing 
the  consistency  of  the  equations,  we  have  retained  it,  since  without  it  the 
equations  would  lack  Gallleian  invariance. 

In  the  Introduction  we  raised  the  question  of  whether  points  on  the  free 
boundary  3fi?f  would  generally  be  regular  points  for  the  Poisson  equation  (1.6). 
This  was  one  of  the  reasons  for  our  search  for  u formulation  of  the  problem 
which  did  not  entail  the  solution  of  a partial  differential  equation  in  the 
liquid  region  subject  to  boundary  data  on  3^.  In  our  reformulation  of  the 
problem,  the  one-phase  Stefan  problem  (3.14)  and  (3.17b)  take  the  place  of  (1.6) 
and  (1.3c).  We  have  indicated  that  (3.14)  is  solved  in  practice  by  using  the 
algorithm  (3.12)  with  da^x')  - da.  The  result  of  such  an  algorithm  has  been 
proven  to  converge  to  a solution  of  the  problem  (3.14)  for  any  given  a * a^  as 
da  -*•  0 (Ref.  2),  and  it  is  not  hard  to  extend  this  to  a proof  of  convergence  to 
the  steady  state  solution.  Our  remarks  in  this  section  indicate  that,  inside 
a boundary  layer  of  thickness  0(T2)  near  the  boundary,  we  may  expect  (2/t2)  v and 
^2/t2)v  to  converge  to  P and  VP,  respectively. 

Accordingly,  it  is  of  some  interest  to  inquire  to  what  extent,  when  the 
boundary  layer  becomes  infinitely  thin  and  the  steady  state  Stefan  problem 
reverts  to  a linear  elliptic  boundary  value  problem,  the  boundary  conditions 
demanded  by  (1.3c)  at  the  free  boundary  are  actually  attained  by  the  result 
of  our  algorithm  in  the  limit  da  0.  Note  that  this* limit  of  an  infinitely 
thin  boundary  layer  can  also  be  achieved  by  letting  Pq  ■*  00  in  (3.11),  (3.14), 
(3.17b),  and  (4.3).  An  error  bound  (Ref.  10)  indicates  that,  for  smooth 
boundaries  3(5^,  the  steady  state  G(da)  given  by  the  algorithm  (3.12)  with 
da1(x')  - da  and  (3.17b)  in  the  limit  ■*•  *>  has  an  L error,  when  N - 3 of 

||G(da)-G(x,x0)||^  - 0 (in  ( )** ),  d - di8t(xQ,  3<?f),  (4.7) 

for  the  computation  of  G(x,Xq)  given  by 

‘o£ 

Gls 


AC-  - 6(x-Xq) 


*0£^f* 


'*f 


- 0. 


<r 
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i 
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(A  recent  paper  (Ref.  3)  in  which  the  same  algorithm  ia  described  appears  to 
give  an  error  O(v^a) , which  is  better  than  the  one  given  above.  I have  not 
tracked  down  the  discrepancy  in  the  estimates.*)  Our  analysis  (Ref.  10)  also 
shows,  for  more  general  regions^,  that  the  steady  state  G(da)  converges  to 
a limit 

C(°)  - G(da)  (4.8) 

as  da  •+  0 , and  that  G(O)-*-  0 as  one  approaches  a regular  point  of  3(f}jU34^)  from 
the  interior  Similar  results  would  apply  to  the  limit  as  da  0 of  the 

steady  state  for  ‘►“computed  by  (3.12)  with  da^(x')  ■ da,  (3.17b),  and  (4.2), 

that  is,  P - 0 at  all  regular  points  of  3«$f  3<$f ) . From  a physical  point  of 

view,  this  may  be  more  reasonable  than  requiring  P ■ 0 at  all  regular  points 
of  3^,  as  in  (1.3c). 

5.  Other  Versions 

We  have  noted  that  there  is  some  arbitrariness  in  the  algorithm  presented 

in  Section  3 regarding  the  presence  of  higher  order  terms  in  T,  as  there  has  to 

be  for  anything  short  of  an  exact  solution  for  that  time  interval.  What  terms 

are  added  or  dropped  is  largely  a matter  of  taste.  For  example,  the  term 

A(uv)  was  left  in  (3.17a)  to  guarantee  Gallleian  invariance. 

Eugene  Isaacson  has  called  our  attention  to  a general  approach  for  solving 

equations  subject  to  a constraint  (Ref.  5).  In  this  approach  one  would  add  a 

perturbation  to  the  equations  so  that  the  constraints  might  be  satisfied,  and 

then  try  to  minimize  the  perturbation  in  some  appropriate  sense.  For  example, 

if  one  had  the  constraint  (1.1b)  on  a velocity  field,  one  could  add  a vector 

field  to  the  right-herd  side  of  the  hyperbolic  conservation  laws  (3.1).  It  is 

o 

well  known  that  the  L minimum  vector  field  with  given  divergence  and  suitable 
homogeneous  boundary  conditions  is  a gradient,  and  one  might  in  this  way  arrive 
at  (1.1a).  In  our  case,  the  only  appropriate  extension  of  this  to  the  steady 
state  one-phase  Stefan  problem  we  are  aware  of  appears  to  be  that,  if  X is  a 

t 

vector  field  satisfying 

V • X - pQ  - p , xe£f , (5.1a) 
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then  the  minimum  of  f^\2 dx  over  all  A satisfying  (5.1)  and  over  allfPjUg  la 
achieved  for 

A « 7v  (5.3) 

+ N 

with  v given  by  (3.17b)  and  (3.14)  with  A replaced  by  A form.  We  have  not 
followed  this  line  of  thinking  in  searching  for  a generalized  hydrodynamics, 
because  we  have  felt  more  secure  proceeding  on  grounds  with  a more  direct  con- 
nection to  physics. 

Another  Item  which  we  have  not  yet  discussed  relates  to  the  degree  to 
which  the  various  stages  of  our  algorithm  conserve  energy.  It  is  obvious  from 
(1.4)  - (3.7)  that  energy  cannot  increase  in  the  first  half  of  our  split-step 
scheme : 

/(*JP u2  + pgz)  d x < / (%pu2  + Pgz)dx  (5.4) 

For  the  second  half,  the  result  is  not  so  simple.  From  (3.17)  and  (4.3)  we 
may  derive 

pu2  - pu2“-p(u-u)2  + ^ v V * u - 2v(Vu)2 

+ A (vu2)  - ^ V • (vu)  . (5.5) 

4 - 

The  only  term  which  can  increase  the  energy  here  is  — v V . u.  To  first  order 
in  T,  we  have  seen  that  an  initially  divergenceless  velocity  field  remains  so 
as  long  as  the  flow  is  a classical  one,  and  hence  to  first  order  there  is  no 
energy  change.  Even  if  we  had  the  possibility  of  a velocity  field  with  a 

<\j 

divergence,  we  would  generally  expect  to  lowest  order  in  T that  p < Pq  where 

V • u >0  and  p > pQ  where  V • u < 0.  Then  to  lowest  order  we  would  expect  v ■ 0 
where  V • u > 0 and  v > 0 where  7 • u < 0.  To  lowest  order  in  x,  we  could  replace 

V • u by  V • u and  thus  conclude  that  ~ v V*  u would  tend  to  reduce  the  energy. 

Our  general  feeling  is  this:  It  may  be  acceptable  for  energy  to  be 

created  in  the  second  half  of  a time  step  as  a compensation  for  too  much  energy 

dissipated  in  the  Inelastic  collisions  of  the  first  half,  but  it  is  not  physi- 
cally acceptable  for  energy  to  be  gained  overall.  On  the  other  hand,  a slight 
energy  increase  which  is  a manifestation  of  the  time  discretization  in  the 
algorithm  as  opposed  to  a sign  of  instability  may  not  be  disastrous.  As  it  is, 
we  can  give  an  example  of  a flow  for  which  the  algorithm  of  Section  3 will  pre- 
dict a net  increase  of  energy  over  a time  step  (Ref.  10).  Needless  to  aay, 

such  a flow  does  not  exhibit  any  great  degree  of  regularity  over  the  time  step, 
and  its  treatment  by  (3.4)  - (3.7)  over  the  first  half  of  the  time  step  Is 
questionable. 


t 


1 We  may  add  a term  to  the  right-hand  side  of  (3.17a)  which  preserves  mo- 

I 

men turn  conservation,  Gallleian  Invariance,  and  consistency  to  lowest  order  in  T 
For  example,  we  may  write 

pGt  - pu±  + A(5tv)  - \ 4.  £5^  Gij  • (5.6) 

In  place  of  (5.5)  we  get 

pu2  - *pu2  - - p(u  -u)2+^-v7*u-2v  (Vu)2 
+ A(vu)2  - -p  V*  (vu) -2  + 2 (UjG^)  • (5.7) 

3 j 

A more  detailed  investigation  (Ref.  10)  suggests  that  is  Is  possible  to  choose 
so  as  to  preserve  Gallleian  Invariance  and  consistency,  and  to  make 

/pu2dx  */pu2dx  1 

in  all  cases.  However,  the  price  may  be  to  replace  (3.17a)  by  a nonlinear 
equation  for  u,  whose  solution  may  have  to  be  obtained  Iteratively.  We  have  J 

gone  to  some  lengths  to  discuss  to  what  extent  energy  conservation  or  nonconser- 
vation Is  an  essential  part  of  our  theory  because  of  our  belief,  elaborated  on  i 

more  fully  In  the  next  lecture,  that  the  status  of  energy  conservation  for  the 
limiting  flow  obtained  as  T 0 has  a deeper  connection  to  Important  qualitative 
properties  of  the  flow  (Ref.  9).  As  we  have  already  indicated,  (3.17a)  is 
likely  subject  to  further  emendation,  and  it  may  be  that  there  is  no  uniquely 
simple  and  acceptable  formula  for  pu,  in  contrast  to  equations  (3.14)  and  (3.15) 
for  P. 

6.  Stratified  Flow  ind  Transonic  Flow 

Incompressible  flows  with  a non-constant  density  are  amenable  to  a treat- 
ment like  that  offered  here  for  the  constant-density  case.  In  this  case  we  1 

Introduce,  in  addition  to  p and  u,  a new  dependent  variable^  which  represents 
the  volume  fraction  of  space  filled  at  each  point.  Our  hyperbolic  conservation 
lews  consist  of 


(/p)t  + V-tfpu)  - 0 , (6.1) 


i 


I 


■ 


f 


As  before,  we  proceed  for  the  first  half  time  step  from  , p,  pu)  to 

(7,  p,  pu)  by  "solving"  (6.1),  (6.2),  and  (6.3)  for  a time  T.  To  satisfy  the 

constraint  (6.4)  we  solve 


with 

Here 


Ba*m  Af*(0*) 

0*(a-O)  - % . 

f*(0*)  - j6*  ' 1 6‘  * 1 

(O  0**1 


(6.5a) 


(6.5b) 


(6.6) 


7 Is  given  by 

f-2£e.<«).  (6.7) 

We  define 

v*  - /“  f*  (0*)da  (6.8) 

and  find  p from 

7P  *^p  + A ( pv* ) (6.9) 

In  place  of  (3.17a)  we  have 

/pu  ■ 7Pu  -^V(pv*)  + A (puv*)  (6.10) 

The  algorithm  in  Section  3 is  then  a special  case  of  (6.1)  - (6.10).  Note  that 
in  this  section  p refers  to  an  intrinsic  flultj,  property , whereas  earlier  in  this 
paper  p refers  to  mass  density.  To  get  the  algorithm  of  Section  3 from  (6.1)  - 
(6.10),  take  the  case  where  p - p - p - pQ  in  (6.1)  - (6.10),  and  then  replace 
7Pq,^Pq,  and^pQ  wherever  they  occur  in  (6.1)  - (6.10)  by  p,  p,  and  p,  respec- 
tively. 

The  equation  of  state  (2.2)  may  be  considered  to  be  a special  case  of  the 
more  general  equation  of  state  for  a barotropic  flow: 


P - P(p).  (6.11) 

Since  we  have  solved  (2.2)  by  solving  a one-phase  Stefan  problem  (3.14),  one  may 
ask  if  (6.11)  can  also  be  obtained  through  the  solution  of  a nonlinear  parabolic 
equation.  One  might  even  wonder  if,  for  a polytropic  fluid,  the  analog  to  the 
one-phase  Stefan  problem  is  the  equation  for  flow  in  a porous  medium.  The 
answer  to  this  seems  quite  clearly  to  be  "No."  We  shall  indicate  some  analogies 


•*-  --  “-“rTTr 
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which  may  carry  over  to  the  transonic  case,  but  we  do  not  think  these  have  any 
practical  computational  value. 


(C)  = sup{x|P(x)  - 5} 


(6.12) 


In  the  hydrodynamic  case,  ^ * (£)  ■ pnV  £ 2 0.  Then  an  analog  to  (3.14)  may 


be  written  as 


eQ  - Af(0(a),  v(o0) 


(6.13a) 


where 


0(oi  - 0)  - p 


f(8,v)  “ max 


v(a)  - £ f(0(a'),  v(a'))da' 


(6.13b) 


(6.14) 


(6.15) 


In  the  cases  of  greatest  Interest, 3^ (0)  • Oand — ^ is  bounded  and  ncm- 
negatlve.  In  such  cases  we  can  show  that  (Ref.  10)  0 and  v achieve  limits  p and 
v as  a that 


Av  ■ p 


2 - 


(6.16) 


and  that  v and  p depend  monotonlcally  on  p.  Also,  If  p has  compact  support,  so 
does  v when  we  have  a polytronlc  fluid 


P(p)  - ApY 


(6.17) 


with  Y > 1* 


Acknowledgement 


This  work  has  been  supported  by  the  Office  of  Naval  Research  under  Task 


No.  NR  334-003. 


1 

1 I 


• ! 


References 

1.  Alan  E.  Berger,  Melvyn  Clment,  and  Joel  C.  W.  Rogers,  The  Alternating  Phase 
Truncation  Method  for  Numerical  Solution  of  a Stefan  Problem,  to  appear  In  SIAM 
J.  Numar.  Anal. 

2.  Haim  Bretls,  Alan  E.  Berger,  and  Joel  C.  W.  Rogers,  A Numerical  Method  for 
Solving  the  problem  u£  - Af(u)  ■ 0,  to  appear. 


- 21  - 


3.  Alexandre  J.  Chorin,  Thomas  J.  R.  Hughes,  Marjorie  F.  McCracken,  and 
Jerrold  E.  Marsden,  Product  Formulas  and  Numerical  Algorithms,  Comm.  Pure  Appl. 
Math.  31,  205  (1978). 

4.  mm—  An  Introduction  to  Probability  Theory  and  Its.  Applications, 

Vol.  II,  John  Wiley  (1966). 

5.  Eugene  Isaacson,  Integration  Schemes  for  Long  Term  Calculation,  Advances 
in  Computer  Methods  for  Partial  Differential  Equations  II,  AICA  (1977) . 

6.  Joel  C.  W.  Rogers,  Water  Waves:  Analytic  Solutions,  Uniqueness  and  Continu- 
ous Dependence  on  the  Data,  Naval  Ordnance  Laboratory  NSWC/WOL/TR  75-43  (1975). 

7.  J.  C.  W.  Rogers,  An  Algorithm  for  a Hyperbolic  Free  Boundary  Problem,  The 
Johns  Hopkins  University  Applied  Physics  Laboratory  APL/JHU  TG  1309  (1977). 

8.  Joel  C.  W.  Rogers,  Steady  State  of  a Non-Linear  Evolutionary  Equation, 
Seminaires  IRIA,  Analyse  et  Controle  de  SystWs  (1978). 

9.  Joel  C.  W.  Rogers,  Stability,  Energy  Conservation,  and  Turbulence  for  Water 
Waves,  Seminaires  IRIA,  Analyse  et  Controle  de  SystWs  (1978). 

10.  Joel  C.  W.  Rogers,  Water  Waves  — Conservation  Laws  and  Constraints,  in 

preparation. 

11.  Joel  C.  W.  Rogers  and  Stanley  Favin,  Computation  of  Water  Waves,  in  prepara- 
tion. 

12.  J.  J.  Stoker,  Water  Waves.  Interscience  (1957). 


